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Plasma flow vs. magnetic feature-tracking speeds in the Sun 
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ABSTRACT 

We simulate the magnetic feature tracking (MFT) speed using advective-diffusive transport 
models in both one and two dimensions. By depositing magnetic bipolar regions at different 
latitudes at the Sun's surface and following their evolution for a prescribed meridional circu- 
lation and magnetic diffusivity profiles, we derive the MFT speed as a function of latitude. We 
find that in a one dimensional surface-transport model the simulated MFT speed at the surface 
is always the same as the meridional flow-speed used as input to the model, but is different in 
a two-dimensional transport model in the meridional (r, 0) plane. The difference depends on 
the value of the magnetic diffusivity and on the radial gradient of the latitudinal velocity. We 
have confirmed our results with two different codes in spherical and Cartesian coordinates. 
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. 1 INTRODUCTION 

At the solar surface magnetic features are observed in the form of 
active regions. Tracking the motion of such a structure individu- 
ally, one finds in general a poleward migration which is suggestive 

■ of a poleward meridional flow at and just beneath the Sun's sur- 
face. However, Doppler measurements of the poleward flow speed 

, at the surface reveal a systematic difference from the speed in- 
ferred from magnetic feature-tracking (MFT): at low and mid lat- 
itudes, the latter is observed to be lower than the Doppler speed , 

I but similar to it at high latitudes (see Fig. 10 of luirichi 20100. 

' Sunspots are usually discarded in such an analysis jKomm et al.l 

, I1993I : iHathawav & R ightmire 2010) because their motion may be 
affected by their strong m agnetic fields. To und erstand the physical 
origin of these differences. lDikDati et al.H2010l) performed a simple 
test using a 2D (axisymmetric) advective-diffusive flux-transport 
model. They showed in simulations that, due to diffusive transport, 
the MFT speed can indeed be different from that of the meridional 
flow fed into the model. They attributed this difference to the lati- 
tudinal gradient of the radial component of the magnetic field, di- 
rected towards the equator at the equatorward side of an active re- 
gion and towards the pole at its poleward side. They concluded that 
magnetic features drift poleward with a net speed that is lower than 
the flow speed at low latitudes and higher at high latitudes. 

In non- axisymmetric two-dimensional (ff, 0) surface-transport 
models (e.g. lSaumann et alj2004ISheelevl201ol ; IWang et al.l2()09l) 
one could likewise suppose that diffusion is the only agent that 
can prevent magnetic features from simply being advected with the 
meridional flow. However, differences between Doppler and MFT 
speeds have never been discussed for those models. 
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In this paper, we use both ID and 2D advective-diffusive flux 
transport models to clarify to what extent the value of the magnetic 
diffusivity and its radial gradient influence the difference between 
the circulation and MFT speeds. Moreover, we will study the role 
of the radial gradient of the latitudinal flow velocity. 



2 MODELS AND METHODS 

For the sake of simplicity we consider the evolution of azimuthally 
averaged, purely poloidal, i.e., meridional fields. Quite generally, 
studying averaged fields requires the inclusion of an additional 
mean electromotive force £ in the induction equation. Its major 
constituents are often described by the a effect, turbu lent pump- 
ing and turbulent diffusivity rjt (see, e.g.. lMoffattlll978i) . However, 
only the latter will be taken into account here and assumed to be 
isotropic, yet possibly depending on depth. We do not claim that 
all other mean-field effects, in particular turbulent pumping, are 
negligible, but prefer to clarify the origin of the speed deviations 
in question by considering the effects in isolation. Thus, we focus 
here on the competition between diffusion and advection. 

Our model is kinematic and we consider axisymmetric solu- 
tions of the induction equation in spherical coordinates (r, 6, (j>) 



dB 

'dt 



V X (t/x B-t/tV X B), V-B = 0, (1) 



with the total (molecular plus turbulent) magnetic diffusivity rjT ~ 
rj + rjt and the prescribed velocity U, i.e., we disregard the back- 
reaction of the magnetic field onto U. The computational domain 
spans over a spherical half-shell _Ri,^r^i?, 0^0^ 7r/2 (i.e., 
from the pole to the equator), where R is the solar radius and the 
base of the convection zone is at Rb = 0.7-R. The total diffusivity 
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Figure 1. Solid lines: Radial profile of Ug{r, 57°) for different k as indi- 
cated in the legend, n = 0.8, cf. Eq. J5). Broken lines: Cartesian velocity 
profile, Uy{x, Ly/2), cf Eq. QT). 



?7t is, unless specified otherwise, constant across the domain and 
considered a free parameter of the model. 

For the ID version of the model, we solve the radial part of 
Eq. l[T} fo r Br- at r = R, as done in several surface-tr ansport models 
(see, e.g.. lDeVore et alJl984l : lBaumann et alj2004l) : 



dBr _ 1 d 



UgBr 



Vt dBr 
R do 



(2) 



Note that this equation is subject to the si mplifying assumptio n 
Be -^ Br at the surface (see the Appendix of lOeVore et al.ll984h , 
the consequences of which will be assessed later when discussing 
our results. We solve Eq. lO by using a second order finite dif- 
ference scheme with 512 grid pointo The time integration is per- 
formed with an implicit (Crank-Nicholson) method. 

For the 2D version we solve instead of Eq. ([T} the correspond- 
ing equation for the (j) component of the vector potential A = Ae^, 
where B = V x A, and 



- = --iU.V)isA) 



■ffr 






A, s = rsinO, (3) 



again utilizing finite differences (Lax-Wendroff scheme for first and 
centered second order scheme for second derivatives). For all sim- 
ulations we use 400'^ grid points. A convergence analysis revealed 
that for the global evolution of the magnetic field a resolution of 
128^ grid points is already sufficient. However, a smoother profile 
of the estimated tracer velocity is obtained with the higher resolu- 
tion. Time integration is d one with the ADI method (for details see 
iGuerrero & Munozll2004h . 

In modeling the meridional velocity U we start with the cor- 
responding mass flow pU which is assumed steady and has thus to 
obey V • {pU) — because of mass conservation. Hence it can be 
derived from a strea m function ip by pU = V x (-i/je^). Following 
iDikpati et alj J2010l) , we assume an adiabatic density profile 



p{r) = poiR/r ^ 0.97y 

where po is specified such that p{R) ~ 5 
stream function we choose the ansatz tp = - 



(4) 

10"^ g/cm^. For the 
tjjoF(r)dgG{e) with 



F{r) 



0.97 



1 



r) 



)'-( 



r \'^ 

r) 



G{e) = P2 (cos 6') + mPi (cos 6) , 



(5) 



(6) 



where Pi is the Legendre polynomial of order I. F{r) guarantees 
vanishing [/,■ at the boundaries r = Rb,R assumed impenetrable. 



Same results are obtained for resolutions from 128 to 1024 grid points. 
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Figure 2. Bipolar region according to Eq. (7) for Q\ = 45°. S± is the po- 
sition of the centre of the polewards (positive) and equatorwards (negative) 
spot, respectively. Color coding: Br . 



The first factor in F{r), resembling the density profile, is neces- 
sary to avoid local extrema of Ue within the domain which can be 
achieved byn = 0.8 ... 1.2, depending on the value of k. Apart 
from that, the exponents n and k are free parameters of the model 
determining the value of the radial derivative of Ue . For fixed n, 
lower (higher) values of k result in a flatter (steeper) radial profile 
Ug{r) (see Fig.[T). Tuning the exponent n allows adjusting the Ue 
gradient just at the surface without changing it very strongly deeper 
down. So n can be employed for ensuring the stress-free boundary 
condition {drUg){R,9) — Ue{R,0)/R, usually imposed in, e.g., 
mean-field hydrodynamic models of stellar rotation and in direct 
numerical simulations of convection. Here, it can be expressed in 
the form hr{R) = 1, where hr{r) is the normalized radial gradi- 
ent of Ue, hr{r) = r{d\nF/dr). For simplicity, we ignored this 
condition in most of our calculations and fixed n — 0.8. However, 
we have checked the influence of having the stress-free condition 
obeyed on our results in a number of cases with different values of 
n. 

Further, the choice of m — —0.2 results in a latitudinal sur- 
face profile Ue (R,d) which resembles Doppler velocity observa- 
tions (see, e.g., IUlrichll2010r ), in particular the position of the sur- 
face maximum of Ue{R, 9) is fairly well reproduced, tpo is adjusted 
such that this maximum is Uo ~ 2500 cm/s. 

The initial magnetic field of a bipolar region is modelled as a 
flux loop in a meridional plane corresponding to two rings of con- 
centrated magnetic flux on the surface of the sphere. We describe it 
by the vector potential 



A = Aq exp 



we 



cxp 



R 



e,/,, (7) 



where 6i is the initial latitudinal location of the center of the bipolar 
region. The initial separation between the positions at which Br 
assumes its extrema at the surface, that is, the "spot separation" is 
<;\/2 We, whereas the depth to which the loop extends is controlled 
by Wr. We assume we — 0.02 (2.3°) and Wr = 0.047? throughout 
this paper. For the corresponding field geometry see Fig.|2] 

Magnetic boundary conditions are chosen to be consistent 
with a perfect conductor in r ^ Rb and to ensure continuity of 
B with an external potential field atr = R. Comparisons with the 
simpler normal field condition r x B = instead of the potential 
field condition showed no noticeable difference in the results. At 
the equator B is assumed to be antisymmetric. 

Both models were run over the model time interval T, being 
typically two weeks, for 20 equidistant initial latitudes 9i of the 
bipolar region between 5° and 85° . For measuring the latitudinal 
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Figure 3. Tracking velocities for the ID model (2) with ijt = lO^^cm^/s. 
Solid line: flow velocity Ug{R), filled and open circles: V+ and V— , re- 
spectively; diamonds: average V = {V+ + V—)/2. 



surface drift velocity of the flux loop, averaged over T, two dif- 
ferent methods were employed. In the first one the position ^q, at 
which the normal magnetic field Br vanishes, was followed. In the 
second, we trace the positions of the local maximum and minimum 
of Br within the loop, 9+ and 9-, respectively. The averaged lat- 
itudinal velocity was then defined as Fo = R{9o{T) - 6lo(0))/r 
in the first case and as the average of the two values V± = 
R {e± (T)~9±{0))/T, that is. 



V = {V+ + V-)/2 



(8) 



in the second. We assign Vo to the average colatitude (^6o{T) + 
6lo(0))/2, but V to the average (61+ (T) + 9-{T) + 9+{0) + 
9- (0)) /4. As the profiles 17(61) and Vq (9), obtained directly in this 
way, turned out to be rather wiggly we fitted them just to Ue{R, 9), 
that is, to the function G in Eq. ^, but with an adjustable ampli- 
tude as fit parameter. For the highest diffusivity used and for start- 
ing latitudes 9i closest to the equator, the influence of this reflecting 
boundary becomes noticeable. This influence leads to an unrealis- 
tically low velocity of the equatorward (negative) spot which, in 
turn, corrupts V. We use instead Vo there. 



3 RESULTS 

3.1 One-dimensional model 

In simulations with high diffusivity, tjt ~ lO^^cm^/s, we find that 
the speed of the poleward spot, V+, is larger than the fluid ve- 
locity, whereas that of the equatorward spot, V-, is smaller; see 
Fig. [5] However, the average velocity V matches the fluid velocity 
Ue fairly well. The velocity of the center of the bipolar region, V^o, 
also coincides with it. For even higher values of tit, V+ and V- 
deviate stronger from Ue, but the average continues to agree with it. 
For evolution times shorter than two weeks (e.g., one week or less) 
the curve for V is more wiggly. However, it always follows the flow. 
These results agree with those of the 2D (9, (f)) model of Wang et al. 
( 120091) , where the poleward spots of the bipolar regions move faster 
than the fluid for a similar value of rfr (see their Fig. 15). 



3.2 Two-dimensional model 

Next, we study the evolution of a two-dimensional bipolar region 
by solving Eq. (O using Eq. (O as initial condition. As a represen- 
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Figure 4. Magnetic field (solid lines) after two weeks of evolution of a 
bipolar region initially at d\ = 45°, see Fig.|2] Color coding: Br- Left: 
??T = lO^cm^/s; right: J7x = 10^'^cm^/s (right). 
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Figure 5. Tracking velocity V(9) for 10*cm^/s <»7t ^3 X 10^-^cm^/s, 
k = 6. Solid/black: flow velocity Ug{R,6). Symbols: V according to 
Eq. (8), omitted for r^-p = 10*cm-^/s as mostly coinciding with those for 
rj^ = 10i"cm2/s; symbol for r?T = 3 X lO^^crnVs, 9 = 85° shows Vo- 
Lines: data fitted to Eq. (6) with amplitude as fit parameter, rra = —0.2 
fixed. Colors/symbols/line styles correspond to different values of r^-p ac- 
cording to legend. 



tative case. Fig. |4] shows the evolution for a bipolar region initially 
located at 9i = 45°, using fc = 6 with either tit ~ lO'cm^/s (left) 
or lO^^cm^s (right). 

In contrast to the ID model, where the profile V{9) turns out 
to be independent of t^t, we find here a significant dependence. 
For small values of t^t, the "frozen-in" condition is fulfilled and 
thus the magnetic field lines appear indeed attached to the plasma 
flow. (The systematic offset between V and Ue for tjt — >■ 0, vis- 
ible in Figs. \5\ and [6l is mainly due to the discretization errors.) 
For larger r^T (<; 3 x lO^^cm^/s), however, the diffusion time 
scale becomes similar to or even smaller than the advection time 
scale, and then there is an increasing departure from the "frozen- 
in" state. This becomes clear in Fig. [5] where V{9) is shown for 
10* ^ ??T ^ 3 x lO'^^cm^/s. In general, the deviation from Ug 
increases everywhere with growing tjt, while for each r]T it adopts 
its largest value at intermediate latitudes of ~ 57° where Ue peaks. 
Note that the simple fit based only on the amplitude, using the func- 
tion G from Eq. l|6j, works remarkably well. 

The top panel of Fig. |6] visualizes the dependence of V on the 
radial variation of Ue, i. e., on the index k in Eq. l|6j. We present 
V at the latitude where Ue peaks (9 ~ 57°) as a function of rjT 
for k varying from 0.5 to 10. To minimize the effect of numerical 
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Figure 6. Upper panel: Maximum tracking velocity V vs. r^x for different 
values of k in Eq. J6). Solid lines/filled circles correspond to the spherical 
model (values taken from fit curves). Broken lines/open squares correspond 
to the Cartesian model (C); dashed/filled squares: density according to 
Eq. (4), dash-dot/open squares: constant density. Lower panel, lines/closed 
circles: fractional speed difference, (Uq{R) — V)/Ug{R), vs. the nor- 
malized radial gradient of Ug, hr(r), averaged over r/R = 0.97. . . 1. 
Open circles: stress-free boundary condition ensured in the profile (3) by 
(n, k) = (0.954, 4) for {K) » 12 and (1.194, 10) for (K) ^ 17.5. 



noise we have taken V from fit curves. For k = 0.5 (yellow line), 

V does not depend on j/t up to lO^^cm^/s. Beyond this value, 

V starts to decrease. For increasing k the curves depart from the 
"frozen-in" domain at decreasing values of 77T being as small as 
« 3 X l(fcva?/s for k = 10. The bottom panel of the same figure 
shows the fractional velocity difference, (Ug{R) — V) /Ue{R), as 
a function of the normalized radial gradient of Ug, h^ (r), averaged 
over the interval r = 0.97i? . . . R, where the major part of the 
magnetic flux is residing. Note that for the highest diffusivity, 3 x 
lO^^cm'^/s, V is reduced by «45% at (hr) ~ 33. Employing ou r 
results for interpreting the data given in Fig. 10 of lUlrichl OOICI) . 
we find that their speed reductions of about 30% do occur in our 
model, either for r^T = 3 x lO^^cm^/s and (hr) ~ 20 or for 
rjT = lO^^cm^/s and (hr) ~ 33. In the bottom panel of Fig.|6] 
some results are shown with the profile (|5} adjusted to the stress- 
free boundary condition by fine-tuning of n. Obviously, there is 
only a slight reduction of the velocity difference in comparison with 
the unadjusted profile. 



3.3 Dependency on r]T {r) 

Having examined the influence of the radial profile of Ug on V, 
one must ask whether also the radial profile of tjt has an ef- 
fect. From an observational point of view this profile is unknown. 
Hence, the profiles so far considered in dynamo models are to 
some extent arbitrary. For instance, IPikpati & Gilmanl ( 1200 ih and 
iGuerrero & de Gouveia Dal Find ( 120070 have used a step function, 
with amplitudes of 10^°cm^/s in the bulk of the convection zone 
and a value of « 10^'^cm^/s for supergranular diffusio n within the 
so-cal led near-surface shear layer. On the other hand, IPipin et al.l 
( 1201 ll) considered a mixing length theory (MLT) estimation of r]t . 
Here we consider both a step and an MLT profile defined 
by the following expressions (see the profiles in the top panel of 
Fig.lIJ: 



step I ^s 

r?T = 7?cz H 



1 + erf 



- OMR 
0.02R 



with J7s — 10^'^cm^/s and rjcz = 10 ^ris, and 



MLT , 

Vt = Vrz + 



+ 



J?cz - ?yrz 

2 


1 + erf 


^r-0.71i?Y 
V 0.02R ) 


r?s - r?cz 
2 


1-erf 


V OMR ) 



(9) 



(10) 



where r^rz = 10 cm /s, 77s = 10 cm /s, and rjcz = 10 rjs. 

We have performed numerical experiments with fc = 6 for 
Ug and a fixed surface value riT{R) ~ lO^^cm^/s. The results, 
displayed in the bottom panel of Fig. [7] do not show marked dif- 
ferences between the three diffusivity profiles considered, and the 
MFT velocity profile is roughly the same. However, the separation 
between poleward and equatorward spots is smaller for the model 
with the step profile. The model with constant t^t exhibits an inter- 
mediate difference whereas in the model with the MLT profile the 
dispersion increases. 



3.4 Comparison with Cartesian geometry 

To assess the influence of the curvature in our spherical model we 
have repeated some of the simulations in a 2D (I/x x Ly) Cartesian 
box with aspect ratio (R — Rt) : R 7r/2 and a simplified circulation 
velocity. 



U = -V X {ipe^ 
P 



with {x, y) corresponding to (r, 9), respectively, p set either con- 
stant or to p{x + i?b) from Eq. (|4j, and i/iq again adjusted to yield 
2500 cm/s for the maximum surface velocity. Instead of the vac- 
uum boundary condition at r = i?, the simpler normal field condi- 
tion e^: X i? = was imposed there. These simulations were per- 
formed with the Pencil CODqj, which uses sixth-order explicit 
finite differences in space and third order accurate time stepping 
method. For these runs we use 512^ grid points. 

The results, here with constant 77T, for both choices of p, are 
represented by broken lines in Fig. |6] The corresponding profiles 
Uy{x, Ly/2) are shown in Fig[T]with the same line styles. Clearly, 
these findings agree with those of the spherical model showing that 
the value of rjT, at which the frozen-in condition starts to fail, de- 
pends crucially on the radial (x) gradient of Uy and varies here by 
about two orders of magnitude. 
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Figure 7. Top: Profiles of »?t used. Bottom: Corresponding tracking veloci- 
ties F for r^x (i?) = lO^^cm^/s, fc = 6. Solid/black: flow velocity Ug{R). 
Filled (open) symbols: V+{V—). 



explicitly on d^ Uy , the surface speed of Bx will clearly be modified 
by the fluid motion deeper down in the sub-surface layers. When 
ignoring By from the beginning, however, as done in PeVore et al.l 
(11984), and many surface transport models afterwards, this influ- 
ence will be lost. In spherical coordinates, the induction equation 
for Br exhibits an analogous dependence on Bg, hence the same 
argument is valid. 

Based upon our results for the difference between flow and 
feature-tracking speeds, one might think of inferring the thickness 
of the layer where the flow is poleward. This value, however, would 
depend on the surfac e diffusivit y and on drllg, both of which are 
poorly known. Hathaway] j2011r) has inferred an extreme flow pat- 
tern with a very shallow poleward flow (~ 35 Mm deep) which 
nevertheless can be brought into agreement with our findings, re- 
quiring a large radial gradient of Ue, that is, {hr(r)} <: 20, cf. 
Fig.|6] On the other hand, surface flux-transport models in two di- 
mensions (6, (j)) which disregard the radial derivatives in Br, are 
probably overestimating the importance of advection in their re- 
sults. 
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4 DISCUSSION AND CONCLUSIONS 

Through ID and 2D advective-diffusive models we have investi- 
gated the differences between the surface meridional flow speed 
obtained from Doppler measurements and that inferred from mag- 
netic feature-tracking. In the one-dimensional simulations, the av- 
erage velocity of the magnetic tracers always coincides with that 
of the flow, independently of the value of t]t. In 2D models, on 
the other hand, flow and feature-tracking velocities may diverge 
at higher diffusivities, for which the "frozen-in" condition does no 
longer hold. Further, the difference between these velocities de- 
pends on the radial gradient of the latitudinal velocity: the steeper 
Ue, the larger the difference. Using a different code we have veri- 
fied that these results apply also in Cartesian geometry. To under- 
stand this dependence we refer to the induction equation (for sim- 
plicity in Cartesian coordinates), taken at the surface x = Lx (or 
r — R) where Ux = dyUx = 0, so 



dBx 
dt 

dBy 
dt 



--Q-iUyBx)^-!!^ 



dx'^ 



+ 



^(t/,B. 



' UxBy) + ?7T 



d^Bx 
d'By 



+ 



dPB, 



(12) 
(13) 



Note that Bx is apparently decoupled from By, but at the price of 
being coupled to its second vertical derivative. In the case of small 
rjT (i.e., diffusion time larger than advection time), the evolution 
of Bx is governed by the first, advective, term in Eq. ( I12l >. In this 
case the magnetic field lines follow the fluid velocity locally (see 
curved magnetic field lines in the left hand panel of Fig.|4}. In the 
case of larger t^t (diffusion time similar to or shorter than advec- 
tion time), the diffusion term in Eq. l ll2t plays a significant role in 
the field evolution. The dependence on dxBx can be eliminated by 
the solenoidal condition dxBx + dyBy = 0, by which the coupling 
between the two equations becomes obvious. Because By depends 
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